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Abstract. We consider a system of Gross-Pitaevskii equations in R 2 mod- 
elling a mixture of two Bose-Einstein condensates with repulsive interaction. 
We aim to study the qualitative behaviour of ground and excited state so- 
lutions. We allow two different harmonic and off-centered trapping poten- 
tials and study the spatial patterns of the solutions within the Thomas- 
Fermi approximation as well as phase segregation phenomena within the large- 
interaction regime. 
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1. Introduction 

The first successful experimental realization of Bose-Einstein condensates for 
atomic gases [4], which goes back to 1995, gave rise to various numerical and 
theoretical investigations on the macroscopic equation ruling these phenomena, 
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that is the Gross-Pitaevskii equation (GPE) 
h 2 

hid t %l> = A^ + V(xU + M\ 2 i(j, xgW 1 , > 0. 

2m 

More recently, in 1997, Bose-Einstein condensation for a mixture of two different 
interacting atomic species with the same mass was firstly realized at JILA [26] 
exhibiting a partial overlap between the wave functions. The vector nature of the 
order parameter gives rise to some intriguing structures and dynamics that are 
absent in the single component case. This, again, stimulated various succeeding 
studies of numerical and theoretical nature. For both single or binary condensates 
we refer the reader to [13,30] and to the references therein. Recently, some efficient 
numerical techniques have been developed to compute ground state solutions of 
GPE [6,9], which can be used to investigate the vector case. On this basis, in 
this paper we deal with the rigorous analysis of the spatial configurations for the 
standing wave solutions (ground and excited states) of the system in M. 2 



(1.1) 



2mi 
2m 



Md t 1p2 = -27^AV>2 + V 2 (x U X 2 )lp2 +tf 2 lfr 2 |V'l| 2 V'2 +M 2 |</' 2 |V 2 



for the unknown tpi : M 2 — * C, i = 1,2, where h denotes the (reduced) Planck 
constant and the coefficients §ij > (defocusing case), with i?i2 = $21, are given 
by the formula [15] 

jYi ■ — |— 777, ■ 

#ij = 2iraij— -, aij = aji, i,j = 1, 2, 

771 j 777 j 

being related to the scattering lengths and mi the atomic masses of the two 
species composing the mixture. Considering the 2D case is not restrictive as there 
are various situations where the full 3D system can be reduced to a 2D system 
with suitably modified coefficients (see e.g. Section 2.2 of [5]). The coefficients 
•da and $12 play the role of repulsive intra-species and inter-species parameters 
respectively. As we will see, when $ 12 is sufficiently large, then some interesting 
overlap and spatial segregation phenomena between the wave densities occur. Con- 
cerning the potentials, we let Vi(xi,X2) = ^(wf-^xi — xn) 2 + lu 2 2 {x 2 — Xi2) 2 ) for 
(xi,X2) in R 2 , where > 0, i,j = 1,2. A typical situation is when the ViS 
have the same center, without loss of generality, the origin. On the other hand, 
there are some physical situations reported in literature, which lead to consider 
off-centered potentials. See, for instance, [32], where the vertical direction in the 
potential is not aligned with the symmetry axis of the trap. Similar equations 
have also arisen as governing equations for electromagnetic pulse propagation in 
"left-handed" materials with Kerr- type nonlinearity [19], in the modified Hubbard 
model in the long- wavelength approximation [22,23], in quadratic nonlinear mate- 
rials with suitable phase matching [18] and in nonlinear optics, for instance in the 
propagation of pulses in a nonlinear optical fiber of bi-modal type due to the pres- 
ence of some birefringence effects generating two pulses with different polarization 
directions [24]. For a wide discussion on nonlinear Schrodingcr systems we refer 
the interested reader to [1-3] and to the references therein. 

Let Ti. be the Hilbert subspace of H^R 2 , C) x iJ^K 2 , C) defined by 

H=\(iP u iP 2 )eH 1 (R 2 ,C)xH 1 {R 2 ,C): f V t (xi, x 2 )\A\ 2 < 00, i = 1, 2 I , 
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which is the natural framework for bound state solutions, endowed with the norm 
\\fa,M\n = T,£r f IW>,| 2 +/ ViixuxiMl*, 

~[ lm i JR 2 Jr2 

and consider the total energy functional E : H — >■ M. associated with 

(1.2) E^ 1 {-,t),^ 2 (-,t)) = Y^E i {i> i (-,t))+$ 12 K 2 / \M;t)\ 2 \M;t)\ 2 , 

l= i J ^ 2 

where, for i = 1, 2, we set 

E l (M;t)) = ^- [ |V^(-,t)| 2 +/ V i (x 1> x 2 )M;t)\ 2 + ^ [ \M;t)\ 4 - 
Mm Jr 2 Jr 2 z Jr 2 

By multiplying the first equation of system (jl.l|) by dtipi and the second by 
dtip2-, taking the real parts, integrating and adding the resulting identities, it is 
readily seen that E is constant on the solutions, namely E(ipi(-,t),tp2(-,t)) = 
E(ipi(-,0),ip2(-,0)) for any t > 0. Also, as for the case of the single equation, 
by multiplying the first equation of (|1.1[) by -0i and the second by ty 2 , taking the 
imaginary parts, integrating and adding the resulting identities, it turns out that 
the total number of particles JVj of the i-th species is time independent (preservation 
of the particle number), namely 



(1.3) / M;t)\ =Ni, t>0, i= 1,2. 
Jr 2 

The ground state (or least energy) solution of Ql.ip is a solution with ansatz 

(1.4) ifi i (x 1 ,x 2 ,t) = e- it t<t> i (xi,x 2 ), (xx,x 2 ) G M 2 , i > 0, i = 1,2, 

where the pair (<fri,4> 2 ) is real valued and minimizes functional (|1.2[) constrained 
to conditions (|1.3p (with ^ in place of ipi). As a consequence, the fas solve the 
nonlinear eigenvalue problem in M 2 



(1.5) 



-2^A0! + 7i(xi,a;a)0i + i9iift 2 |^i| 2 0i + i9i2^ 2 |0 2 | 2 0i = Mi0i 

■^A(j>2 + V 2 (xx,X 2 )<f> 2 + l9 2 lft 2 |0l| 2 02 + ^22^ 2 |0 2 | 2 02 = M202 

1 2 



Testing the first equation of (|1.5[) by ^£>i and the second by ip 2 , we have a formula 
for the eigenvalues ^ (also known as chemical potentials) versus the eigenvectors 

(1.6) Nm^EM + ^h 2 f |4| 4 + « 2 / |0i| 2 |0 2 | 2 , i = l,2. 

^ Jr 2 Jr 2 

The existence of ground state solutions to (|1.1|) in 7i is reduced to the existence of 
minima for the energy functional (|1.2|) constrained to 

(1.7) S = Wii 2 )£«: ||0i||i2 = iVi, 2 = 1,2}. 



As the are positive, 



^|0l| 4 +^12|^l| 2 |^2| 2 + ^|02| 4 >O, 
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so that the energy functional E is coercive, bounded from below and weakly lower 
semi-continuous over S. Hence, the existence of a ground state solution is im- 
mediately guaranteed. Usually, with reference to the solutions of the form (|1.4p 
(standing waves), there are two possible (physically different) approaches depend- 
ing on whether one considers the chemical potentials fii as fixed (hence searching 
for solutions to the first two equations in (|1.5|) but with possibly different L 2 norms) 
or the total masses J R2 \4>i\ 2 as fixed (thus solving the nonlinear eigenvalue problem 
(|1.5j) . which is the case we deal with). Any other solution (<fii, (f> 2 ) of system (|1.5j) 
of the form (|1.4p not having minimal energy for E will be called excited state (or 
higher energy solution). 

The main goal of this paper is to prove some geometrical properties (clearly 
confirmed by some numerical simulations) for ground and excited state of (| 1 . 5[) . 
particularly under the influence of strong interaction effects (namely §i 2 — > 00). 
See e.g. Propositions 13.11 and 13.21 

In Section [2] we derive the location of ground state solutions via the Thomas- 
Fermi approximation and classify the relative configuration of 4>i with respect to 
4>j . In Section [3] we study the phase separation process (spatial segregation) in the 
large competition regime by means of suitable limiting energy levels which provide 
K-indepcndcnt upper bounds for the energy of solutions. In ScctionJH for the sake of 
completeness, we briefly describe the functional framework of the numerical scheme 
used to compute the solutions. 

In the following the Hubert space L 2 (R 2 , C) is endowed with the standard scalar 
product (/, 3)2 = J R 2 f§, f,9 € L 2 (R 2 ) and the induced norm is denoted by || • \\ L 2. 

2. Location and Thomas Fermi approximation 

If the distance between the centers of the trapping potentials Vi is sufficiently 
small compared with the radii of the supports of the ground state solutions </>i, 
then the condensates share a region where they coexist (with a partial or full over- 
lap, that is one condensate is partially or entirely included in the other). In the 
opposite case the supports of the wave functions are disjoint. Hence, we can en- 
counter three different patterns for the spatial wave functions (f>i, which we are 
going to discuss, namely: no overlap, partial overlap and full overlap. It should 
be noted that support just means here the planar region where the mass of the 
ground state solution is mainly concentrated, being (exponentially) vanishing on 
the outside. In the Thomas-Fermi regime, an approximation of the ground state 
solutions of system (|1.5| . which is very good for sufficiently large values of the 
coupling constants, can be obtained by simply dropping the diffusion terms — A0i, 
namely the kinetic contributions, thus assuming the wave functions to be slowly 
varying (cf. [16,20,21,35]). In turn, (jl.5|) reduces to the algebraic system (here 
we let h = 1) 

f 2tf n |</>i| 2 + 2i9 12 |0 2 | 2 = 2 Ml - ( Xl - Xll ) 2 - (x 2 - x 12 ) 2 , 
(2-1) ^ 

[2021^1 1 2 + 2$ 22 \q> 2 \ 2 = 2^2 - (Xl - X2l) 2 - (X 2 ~ X 22 ) 2 , 

where the /i;S should be computed through the normalization conditions (|1.3|) (if, 
for instance, d\ 2 = 0, it holds \ii cx v^il f° r i = 1, 2). In general, as the left-hand 
sides are positive, this system is satisfied only on a (possibly empty) subset O CM 2 
(O = 0i n 02 in the notations introduced below), namely the overlap region. It is 
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natural to introduce the circumferences defined by (x±— xn) 2 + (x2 — Xi2) 2 



with 



ri(ni) = \Z2/ii, i = 1,2. The intersection of the corresponding disks X>j gives the 
region where system (|2.1[) makes sense. Outside the region O, the wave functions 
take the usual form of the solutions of the GPE ($12 = 0) in the Thomas Fermi 
regime (see the expressions below of <f>i on \ O). More precisely, we have the 
following non-smooth approximations of the ground state (see also the work by 
Riboli and Modugno [32]) 



2(tfll#22-#? 2 ) : 



2j£i -(x 1 -x 11 ) 2 -(x 2 -x 12 ) 2 



0, 



1)2 = < 



$11 (2/J2 — (xi -3:21) 2 - (X2 -x 22 ) 2 ) — #12 (2jJi -(zi -an) 2 - (^2-a:i2) 2 ) 

2(tfll1?22-lJ 2 2 ) ! 



2^2 -(3:1—3:21) 2 -(^2-3:22) 2 
21J22 



0. 



in O, 

in Pi \ O, 

in M 2 \X>i, 

in O, 

in P 2 \ O, 

inM 2 \P 2 , 



that is, equivalcntly, 



(2.2) 



\/ai(Rf - (xi - jm) 2 - (x 2 - 2fe) 2 ), in 0, 

/ r 2 -(x 1 -Xii) 2 -(x2-x i2 ) 2 • ^ v ro 

V 2?- > mVi\U, 

0, in R 2 \ T>i 



where, according to the notations introduced below, O = 0\ H O2, T> = T> \ n T>2, 
= {(xi,x 2 ) € V : {xx - Vll f + (x 2 - Uaf < (>)-R 2 }, 
Vi = {(xi,x 2 ) G R 2 : (xi - x 4i ) 2 + (x 2 - x l2 ) 2 < r 2 }, 
with < (resp. >) in the definition of Oi for a% > (resp. a% < 0), if we set 



^22^11-^12X21 CJi2 W 2 2^12 - ^12^22 

1/11 = = £11 H Aix, j/12 = 

UJ22 — ^>12 OL\ W22 — W12 

CJ11X21 - W12X11 UJ12 U)uX 2 2 - ^>\2X\2 

2/21 = = ^21 AiX, 2/22 = 

— w\% a 2 lu 11 — lj 12 

where AjX = x\j — X2j, namely, for i, j = 1,2 with i ^= j, 

Vij — Xij 



, Wl2 A 

X12 H A 2 x, 

X22 - — A 2 x, 

OL2 



with LOi.\ = 



0, 



y 2dct e 



s, i, j = 1, 2, ai = w 2 2 - W12, a 2 = W11 - W12, and 



2 2w 2 2Ml - 2^12^2 + ^12X 2 i + ^12X^2 - W 22 X 2 i - W 2 2X 2 2 , 2 , 2 

H i = 1" 2/n + 2/i2 

W22 — &12 

2 2^11^2 - 2^12^1 +UJ 12 X 2 11 +U) 12 x\ 2 - Wlixlx - ^11X^2 l2 ,2 

^2 = 1" 2/22 + 2/21" 

07n — 07i2 

Setting Xi = (xn, xq) and y, = (yn, 2/42) for i = 1,2, this reads as 



rf + ——(^i - Mi) + ^(IxjI 2 - |yi| 2 ) - ^(| Xi | 2 - |y ; | 2 ). 

ai ai ai 



(i 
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Hence, we have four circumferences that rule the geometry of the ground states 

E[ : (xi - xn) 2 + (x 2 - x i2 ) 2 = r 2 , Ef : (x\ - yn) 2 + (x 2 - ?to) 2 = R\- 

If ?9i2 = (namely no interaction), we deduce that Xj = yi, Ri = ri, and Oi = T>i 
($12 = implies on > 0), so that the ground state solutions turn into the usual 
Thomas-Fermi representation for the single GPE 



2&M 



0. 



If $12 « 0, then uoii 
for i, j = 1,2 
yi = Xi but different radii R 2 



l 

2l?22 

so that Ef 



« EJ for i = 1,2. 



x 2 , the Ef-s have centers 



2^12 



(Mi 



0, Q!j 
If Xl 

/ij), for any i 7^ j. 



2.1. Nonoverlap case. In the case occurring when the constant $12 is zero 
(absence of interaction in the mixture), the system uncouples into a pair of GPEs 
(for the single GPE various accurate and efficient numerical techniques have been 
recently compared in [9]). Numerical experiments show that the ground state 
solution <pi always locates its mass around the minimum point Xi = {xn, x^) of Vi, 
for i = 1,2. It looks apparent that boosting up the parameter da in front of the 
cubic nonlincarity in the equation of <pi has the effect of squeezing down the profile 
of <pi making it flatter and larger. Going back to the case $12 ^ 0, we say that 
we have no overlap between <f>\ and <p2, if the centers Xi and x 2 of E[ satisfy the 
geometric condition 



(2.3) 



(Ax*) 2 + (A 2 x) 5 



namely if Xi is sufficiently far from x 2 with respect to the amplitudes of the 
supports of 4>i. In this situation the ground state solutions look like those of the 
decoupled case. In fact, the coupling terms CV,- = l di2\4>j\ 2 4'i with i 7^ j arc almost 
everywhere zero as the supports are disjoint, due to (|2 .3[) . Hence the system is 
actually a small deformation of a pair of uncoupled GPEs. Sec Figure [TJ 





10 20 30 40 



Figure 1. Supports of <\>\ and (f>2 (starred) and disks (unstarred) 
related to the Thomas-Fermi approximation for 4>i and (f>2- We 
have taken Ni = N 2 = mi = m 2 = 1, $11 = 400, $22 = 200, $12 = 
100, xn — 4.5, x 2 i = —4.5 (left figure, tangential supports); and 
xn = 6, x 2 i = —6 (right figure, disjoint supports), with x^ = 
for any other 
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2.2. Partial overlap case. We have partial overlap between fa and fa, if 

namely the disks of boundaries E^ and E!J overlap without being one completely 
embedded in the other. We see in the contour plot of Figurc[5]thc partial overlap in 
the case $n 3> $22- Apparently, boosting i?n with respect to $22 makes the overlap 
region more localized. Keeping in mind the behaviour of the uncoupled case, in 
order to give this fact a very simple empirical explanation, it suffices to argue on 
the coupling terms CV,-. In the region (depending upon the relative magnitude of 
the i?jjs) where both fa are nonzero the contribution of Cij pushes down the profile 
around the origin (the center of trapping for fa ), provided that $12 is significantly 
large. The support of fa still remains contractible, but the radial symmetry prop- 
erty of fa is broken (due to strong interaction). See e.g. the situations reported 
in Figures [5] and [5l As Figure 2] shows, while the Thomas-Fermi approximation 
disks E^ are overlapped to the support disks E[ when the coupling constant $12 is 
much smaller than the $us, in the case where $12 3> i.e. in the large interaction 
regime, the four circumferences intersect in two points and E[ and E^ may have 
quite different sizes. 

8 

6 

4 

2 


-2 
-4 
-6 

'ho -5 5 10 15 

Figure 2. Supports of fa and fa (starred) and disks (unstarred) 
related to the Thomas-Fermi approximation for fa and fa (partial 
overlap case). We have N\ = N 2 = m\ = m 2 = 1, i?n = 400, 
^22 = 200, $12 = 100, x\\ = 2, X21 = —2 and Xij = for other 

2.3. Full overlap case. We have full overlap between fa and fa, if 

(A^) 2 + (A 2 x-) 2 < 1^277- T^l 2 , 

so that the disks of boundaries EJ and Ej are included one in the other, see Figure[3] 
In a highly interacting regime, this configuration leads to a non- contractible support 
for one of the two wave functions, which looses the symmetry properties of the trap. 
If the potentials are both centered at the origin, # n ^> #22 > and the coupling $ 12 is 
sufficiently large, as fa spikes around the origin, fa feels the influence of the coupling 
^12 1 03 i\ 2 fa, lowing down the profile (around the origin) and giving rise to a local 
minimum. This behaviour will be rigorously justified in the forthcoming section via 
energy estimates. In the Thomas-Fermi regime the location of the overlap regions 
depends upon the values of dij and of the centers Xi according to formula (|2.2jl . In 
order to visualize some situations arising with respect to the position of the trap 




8 



MARCO CALIARI AND MARCO SQUASSINA 



centers, see Figures [6] and [7] where we kept the $ijS fixed and varied the position of 
xi (xn = ±2 in Figured] and %12 = ±2 in Figure[7|), while X2 = 0. 




Figure 3. Supports of <fii and 4>2 (inside of starred disks) and disks 
related to the Thomas-Fermi approximation for <p\ and (f>2 (overlap 
inside the smaller disk). We have taken Ni = N2 = mi = to 2 = 1, 
i?n = 400, $22 = 200, ??i2 = 100 and XlJ = for other ij = 1, 2. 




Figure 4. Supports of <f)\ and 4>2 (starred disks) and disks related 
to the Thomas-Fermi approximation for cj>\ and 4>2 (which holds 
within the intersection of the unstarred disks). We have taken 

Wl = 7V 2 = mi =7712 = 1, tf n = 400, ??22 = 200, #12 = 1 

(left figure) and $12 = 150 (right figure). In both figures we have 
x\\ = 2, X21 = —2 and Sy = for other 



3. Strong interaction and phase separation 

In the next sections we deal with the justification of the phase separation phe- 
nomena occurring when the repulsive interaction between the condensates gets 
very strong. We consider both ground and excited state solutions. See also 
[10-12,14,25,28,33] for various studies of spatial segregation phenomena in sys- 
tems with large interactions. 
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3.1. Ground state solutions. Assume that the intra-species parameters fins 
are chosen within a bounded range of values and, on the contrary, that the inter- 
species interaction rate $12 becomes very large, say $12 — K, where we let the 
parameter k > go to infinity. For notational simplicity, we set h = 1. Let 
TL C H 1 ^ 2 ) x iJ^R 2 ) be the realization of the Hilbert subspacc given in the 
introduction and consider the energy functional E K : TL — > ffi re-written as 



(3.1) E K (<t> 1 ,<j> 2 ) = E 0O (ct) 1 ,<j>2)+K 

Jm 2 

where 

2 

Eoc(<t>i,fa) =^E i {<j> i ). 

i=l 

Recalling (JTTTJ) , the energy level of the ground state solutions is 

c K = inf E K ((f>i,(f> 2 ). 

(01,02)65 

We also define the value for the limiting segregated least energy value c OCll 

Coo = inf -Eoo(0i,02), 

(01,02)65=0. 

where we have set 

Soc = {(4>i, 4>i) G S : fafc = a.c. in R 2 }. 

Proposition 3.1. The sequence (</>£ , </)%) c S of ground state solutions of (jl.5|) 
converges in TL to a function (0i°,<?!>2 O ) at the energy level c^. Furthermore, 

(3-2) - ^rA^ + Vi{xx,x^T + M<t>T\ 2 <t>r < l*?4>?, 



wh 



ere 

1 JR 2 

fori = 1,2. 

Proof. The infimum that defines the value Coo is taken over a smaller set with 
respect to the one defining c K . Moreover, if the functions 4>i an d 02 have disjoint 
supports, E K {4>\,4>2) = Eoo{4>\, fa), for any k > 0. In particular, of course, this 
implies that c K < Coo, for all k > 0. Therefore, for the ground state solutions 
(<t>1, G TL, §1 ^ for i = 1, 2, we have E K </>£) = c K and 

xk|2| jki2 ^ 771 /J, re J,re\ 1 ,.~ / I j, re 1 2 1 j re 1 2 



(3.3) «/ iffi^ + « / Wri^r = c « < Co 

JR 2 

for every k > 0. As a consequence, 



(3.4) lim / |ffl^=fl. 

JR 2 

Also, for all k > 0, we have 

||(^,^)||^ < EM,<^) + K [ \<ft\ 2 \<%\ 2 < Coo, 

JR 2 

so that the sequences (4 > ti4>2) IS uniformly bounded in TL. In particular, up to a 
subsequence, there exist (</>i°,<^2°) m W such that (^1,^2) ~^ i'PT^T) m ^ as 
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k — ► oo and (j)f(xi,X2) — > (/>i°(xi, X2) a.e. in R 2 . Hence, by combining Fatou's 
Lemma with formula (|3.4j) , we get 

f i^ri 2 i^i 2 = o, 

Js. 2 

so that 

(3.5) ^fW = 0, a.e. in R 2 . 

Since by definition of ground state solution J R2 \cf>f | 2 = Ni, for any k > 0, and H 
in compactly embedded into L 2 (R 2 ) x L 2 (R 2 ) (see inequality Q3.8|l ). up to passing 
to a further subsequence, we conclude that 

(3-6) / \<t>?\ 2 = Ni, 

Jm 2 

for i = 1,2. In particular (</>i°, </>2°) € Soc, by virtue of (|3.5|) and (|3.6[) . Observe 
also that, by virtue of (|1.6|) . (|3.3|) and Ei(<pf) < c^, 

sup^ = ^sup(i? J :(^) + ^ / |tf| 4 + «/ |^| 2 |^| 2 1 < 00, 

being /i" the eigenvalues corresponding to cj)f. Then, up to a subsequence, fif — > 
as k — > 00. By testing the equations of (|1.5[) by an arbitrary positive function 77 
with compact support, we get 

^- I V0f-Vr?+/ ViixuXiWv + tu f m 2 tfri<ii!? [ tfv, 

for all re > and any 77 e C£°(R 2 ) with 77 > 0. Hence, letting re — » 00, it turns out 
that </>°° satisfies the variational inequality (|3.2|) . Notice that, since (</>i°, </>2°) £ ^00, 
by the definition of Coo , we deduce 

2 22 

^ 2m, 7 R 2 ^ 7 R 2 ^2 Jjj2 K^OO ,/ R2 

2 2 2 

^ E ^- lim inf / i v ^i 2 + E lim inf / ^#"i 2 +E% lim inf / k K 1 4 

+ lim K / \<f>1\ 2 \<t>2\ 2 < lim inf E K = lim inf c K < Cco < , <f>?) 

=E^r/ iwn 2 + E/ «°i 2 +E^/ i^°i 4 > 

^ 2TOi 7 R2 ^ 7 R2 ^2 7 R2 

which yields re f R2 |0i| 2 |02| 2 ~* as K ~* 00 7 which is a much stronger conclusion 
compared with ()3.4j) . Consequently, the convergence of <pf to 0°° in H is strong, 
otherwise, assuming by contradiction that for some i = 1, 2 

/ |V<AH 2 < I™ / |V^| 2 or / K^ff < 1™ / ^|^| 2 , 

the previous inequalities we would become strict, yielding immediately a contra- 
diction. Finally, as a further consequence, c K — -> Coo as re — > 00 and the value is 
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indeed assumed and 




Finally, the strong convergence and (|1.6[) yield 

Z JR 2 

for any i = 1,2, which concludes the proof. □ 

As one can see in the numerical simulations, as the interaction coefficient gets 
large, the phase separation becomes rather evident. Sec Figures [5] and [6][7] (just </>i 
component) where different choices of the centers of the T^s have been considered. 



4>i 4>2 




-10 -5 5 10 -10 -5 5 10 



x x 

Figure 5. 2D contour plots, in the square [—11, ll] 2 , of the 
ground state solution for N\ = N2 — mi = m-2 = 1, $11 = 850, 
"$22 = 18, $12 = 210, where the potentials have centers x\\ = 4 
and Xij = for any other i,j. The phase separation is evident 
around the origin (partial overlap case). 

3.2. The anisotropic case. Depending on the relative magnitude of param- 
eters Wjj, there are some directions along which the ground state solutions tends 
to concentrate. For instance, for uj\\ (resp. OJ12) much larger than uj\2 (resp. uju), 
the component 0i has a cigar- like shape along the y-axis (resp. x-axis). Similar 
behaviour for <p2 along the y-axis (resp. x-axis) for W21 (resp. U22) much larger 
than u>22 (resp. t^2i)- In Figure [8] we consider the small interaction case, namely 
$12 <C when to a = 100 and Uij = 1 for i 7^ j. As it is evident from Figure [9l 
increasing the inter-specific coupling constant (1^12 = 1200) the wave functions </>i 
and 4>2 spatially segregate around the origin. 

3.3. Excited state solutions. As for the ground state solutions, in the strong 
interaction regime, also higher energy solutions exhibit a phase separation be- 
haviour. These phenomena have also been confirmed by some numerical simu- 
lations, see e.g. the comparison in Figure [TOl (top. i?i2 = and bottom, i?i2 = 120) 
and in Figure QT] (top, $12 = and bottom, i?i2 = 120). See also Section |4] for the 
notations. 
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x x 

FIGURE 6. 2D contour plot, in the square [—11, ll] 2 , of the first 
component of ground state solution for N\ = N2 = mi = = 1, 
$11 = 850, i?22 = 18, #12 = 210 in the cases where the potentials 
have centers x\\ = 2 and Xij = (left) and x\\ = —2 and x^j = 
(right). The symmetry breaking is evident (full overlap case). 




-10 -5 5 10 -10 -5 5 10 



x x 

Figure 7. 2D contour plot, in the square [-11, 11] 2 , of the first 
component of ground state solution for N\ = N2 — mi = to 2 = 1, 
= 850, $22 = 18, $12 = 210 in the cases where the potentials 
have centers X12 — 2 and Xij = (left) and X12 = ~2 and Xij = 
(right). The symmetry breaking is evident (full overlap case). 



Consider the energy functional (|3.ip defined on the space TC. If we consider 
the family of all subsets A C H \ {(0,0)} which are closed and symmetric w.r.t. 
the origin, the Krasnosclskii genus of A ^ 0, denoted by j(A) £ N, is defined as 
the smallest positive integer n such that there exists an odd continuous function 
£ : A — > M™ \ {0}. We also set 7(0) = 0. When such an integer n fails to exist, we 
put j(A) = 00. Given a positive integer m, we can now introduce the families £, 
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4>i 4>2 




Figure 8. 2D contour plots, in the square [— 5, 5] 2 , of the ground 
state solution for N\ = N2 = mi = m-2 = 1, $11 = 400, $22 = 150, 
$12 = 1, tun = 100, UJ22 = 100, lu\2 = W21 = 1 and potentials 
centered at the origin. 



01 02 




Figure 9. 2D contour plots, in the square [—5, 5] 2 , of the ground 
state solution for N\ = N2 = mi = m% = 1, $n = 400, $22 = 150, 
i?i2 = 1200, u>n = 100, cj 2 2 = 100, U12 = W21 = 1 and potentials 
centered at the origin. The segregation around the origin is evident. 



£0, r m and r™ of subsets of H, defined as follows: 

£ = {A cH \ {(0, 0)} : A is closed and symmetric w.r.t. the origin}; 

£ = {A e £ : if (0i, 02 ) £ -4 then 0i0 2 = a.e. in M 2 }; 

T m = {Ae£: 1 {A) > m and J R2 2 = N u J R2 0^ = N 2 for any (0i, 2 ) G A); 
r™ = {A € £0 ■ l(A) > m and J R2 2 = 7V 1; / R2 2 = AT 2 for any (0 X , 2 ) G A}. 
Then the candidates values to detect some critical (higher) levels of E K are 

c™ = inf sup S K (0i,0 2 )- 

We also introduce the following K-independent values (notice that E K \e = i?oo), 

c™ = inf sup -Boo(0i,02)- 
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Since So C £ , we have T™ C T m for any m and, then, by the above definitions, 

(3.7) C<c™> for all m e N and k > 0. 

As we prove, the levels c™ (which satisfy c™ < c" i+1 as r m+1 C F Tn for any m G N) 
correspond to critical points of E K on W constrained to the sphere 5, thus yielding 
a sequence of nonlinear excited states of the Gross-Pitaevskii system (jl.5p . 

Proposition 3.2. Let m a positive integer. Then, there exists a sequence of 
solutions (</>"'"\ <^2' m ) of system (|1.5|) at energy levels c™ smc/i i/iat, m i/ie Zarge 
competition limit k — > oo, it converges, weakly in Ti. and strongly in L q (M. 2 ) for any 
q > 2 to a limit configuration (</>J°' m , </>^°' m ) G iSoo. 

Remark 3.3. Contrary to the case of ground state solutions it seems not possi- 
ble to show that the limiting configuration {iff' m , 4>^' m ) corresponds to the energy 
level c™ for the functional Eqq and satisfies suitable variational inequalities. 

In order to prove Proposition 13. 2[ we first show that, since Vi — > oo for 
(xi,X2) — ► oo, E K satisfies a technical compactness condition, the Palais-Smalc 
condition. For the sake of completeness, we shall include a proof of this fact. 

Lemma 3.4. For any k > the functional E K \$ satisfies the Palais-Smale 
condition, namely for any sequence (</>*, </> 2 ) in S such that 2£ K (<^* , <?5^) is bounded 
and dE K \s{(j)ln <Pn) ~ y as n ^ oo in the dual space Ti* of Ti (called Palais-Smale 
sequence) there exists a strongly convergent subsequence in Ti. 

Proof. Let n > and let (</>*, </> 2 ) C S be a Palais-Smale sequence for E K . In 
particular, 

sup||(^,^)||« < supE^l,^) < oo 

n>\ n>l 

Hence (0* , </>^) is bounded in Ti and, up to a subsequence, it converges weakly in Ti, 
and for a.e. (x±, X2) in M 2 , to a function (0^-,, f^) € 7i. Notice that Ti is compactly 
embedded into L 2 (R 2 ) x L 2 (R 2 ) as, for any £ = 1,2, we have 

(3.8) sup sup R 2 / (0,\) 2 <oo. 

Then, up to a further subsequence, 0* — » 0^ in L 2 (R 2 ) as n — > 00, which yields 
(0£o,</>^o) € 5. Hence, by the Gagliardo-Nircnbcrg interpolation inequality 

(3.9) U\\ LT ^ m < c!|V0||£ 2(K2) ||0||^ ( Q R2) , Va e [0, 1), V</> e ^(R 2 ), 
taking, in particular, a = 1/2 we get (c > changes from inequality to inequality) 

110; - < c||V< - H^jll^ - O^V) < - #JIl2( R2) 

so that <^ converges to </>^ strongly in L 4 (R 2 ) as n — > 00 (actually in any L 9 ), for 
i = 1,2. Now, by virtue of the condition d-E^ls^* , 2 ) — > as n — > 00, there exists 
a sequence (if„) in Ti* with io„ — > in Ti* as 77 — > 00 and two sequences (fj, l n ) C R, 
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i = 1,2, such that, for all (ip, 77) € Ti 
1 



2m 1 
1 



V^-Vt?+/ V 2 (x 1 ,x 2 )^ n r 1 + K \K\ 2 <V + #22 IKIKV 



2m 2 

(3.10) = /i* / r \v3 + /4 / 0n77 + (<Pi 7 ?)) 



Observe that, by choosing tp = (fi^ and 77 = (resp. y> = and 77 = <p\) and 
recalling that / R2 (0n) 2 = N i: we get a representation formula for /j* (resp. 7J 2 ). It 
follows that {n % n ) is bounded in R so that, up to a subsequence, it converges to some 
positive number 71^. Finally, choosing any arbitrary (ip, 0) £ Ti and (0,77) G 7i as 
test functions in the previous identity and taking the limit as n — > 00, it holds 

(3.11) 7^-1 v^-v^+/ v 1 (x 1 ,x 2 )<g p + d ll [ \4>U 2 4Lv 



2ra\ 



(3.12) ^- / V<f> 2 QO -V V + [ V 2 (x 1 ,x 2 )<f> 2 00 ri + K [ IfaUloV 

^m 2 J R 2 J & 2 J K 2 

+ 7? 22 / |^| 2 ^ = ML / ^ 

In particular (0^,, € Ti is a weak solution of 

-5^A^ + ^(xi.xa)^ + ^iil^Ll 2 ^ + ^I0 2 ool 2 ^ = /4^c 

-S^A^ + ^(xi,a: 2 )02o + >#Ll Voo + WLlVoo = 

Now, choosing <p = and 77 = ^ in (|3.10[) . (/? = 0^ in (|3. and 77 = (j) 2 ^ in 
()3.12l) . taking into account the strong convergence of <f> l n to c/f^ in L 4 (M. 2 ), that 
(</>*, 0^) is bounded in Ti and w n — > in 7Y*, by the resulting identities we get 

2 

lim ||(0»)ll« = l™ f l V ^| 2 + / ^(^i^2)(C) 2 

rwoo n-»oo ^ Imi J R 2 J R 2 



i=l 
2 



lim 

n — >oc 



N^l + N^l-y^Vu ICI 4 -2 K / |#| a |#| a + <«>«> 

2 

= N 1 ^ 00 +N 2 ^ 00 -^ H f W 00 \ A -2nf |0L| 2 |0 2 oo| 2 

2 

= E^r/ IWoo| 2 +/ ^(*i,s a )(^) a HI(#c,0*. 

JZJ ^ TO « JR 2 JR 2 

where we used the fact that — > (fix^oo m £ 2 (R 2 ), following by 

^ I^n0 2 -^L^Ll 2 < 2 II 1 1 i.-iCIRS^, 1 1 0^. — 1 1 i,J <:n£ 2 ) 2 1 1 || ^-l^sj || <^>^. — lll,-i CIR2) - 

Hence </> 2 ) converges in Ti, proving the Palais-Smale condition. □ 
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We now recall the following existence result (see e.g. [34, Theorem 5.7]). 

Let (X, || • ||) be a infinite dimensional Banach space and let Y C X \ {0} be a 
complete symmetric C 1,1 -manifold. Let / : Y — > R be an even functional of class 
C 1 . Assume that / satisfies the Palais-Smale condition and is bounded from below 
on Y. Then / admits at least N — sup{7(A') : K C Y compact and symmetric} 
critical points. 

We are now ready to prove Proposition 13. 21 

Proof of Proposition \3.Sl Since E K is a C 1 functional, satisfies the Palais-Smale 
condition by Lemma f3 .41 is even and bounded from below (as E K > 0), the above 
mentioned result applies with Y = S yielding (it holds N = oo) a sequence of 
solutions (4>i' m , </>%' m ) in S to (|1.5j) with E K (4>1' m , (f>2' m ) — c™, m > 1 and k > 0. 
With reference to ()3.1[) . by means of ()3.7() this implies that 

/ i ;/t,m |2 1 jK,m |2 ^ 771 / iK.Tn iK.Tn\ , / 1 1 K.m |2 1 iK,m |2 m ^ m 
K / 101 I 102 I <#°o(01 ;02 )+ K / 101 I 102 I = C K < C oo, 
il 2 il 2 

for every k > 0. As a consequence, as c™ is independent of k, for any m > 1, 



(3.13) 




Similarly, as || (<fii' m , 02' m )||« < c™ < c™, it follows that (4>i' m , 02' m )«;>o is bounded 
in 7i. Hence, up to a subsequence, (0i' m , 02' m )re>o weakly converges in 7i (and 
strongly in L 9 (R 2 ) for any q > 2, by combining formulas (|3.8p - (|3.9[) ) to a function 
^oo,rn^oo,m-j ^ ^ j n particular, 

/ \^ m \ 2 =N t and <^< m ^°< m =0 a.c. in R 2 , 
Js, 2 

which proves the assertion. □ 



4. Numerical computation of solutions 

We describe the numerical algorithm used for the computation of the ground 
states for the single one-dimensional Gross-Pitacvskii equation and we mention at 
the end of this section how the same technique can be applied to a system of any 
number of coupled equations in R 2 . Moreover, without loss of generality, we reduce 
to the case H = m = 1. The main idea is to directly minimize the energy E((f>) 
associated to a wave function ip( x ) — e~ v '0(x), discretized by Hermitc functions. 
As it is known, the Hermite functions (Tif)i^ are defined by 

H?(x) = E^{x)e-^ x \ I E N, 

where (H^)i^ are the Hermite polynomials [7], orthonormal in L? with respect 
to the weight e~@ x . The Hermite functions are the solutions (ground state, for 
/ = 0, and excited states, if else) to the eigenvalue problem for the linear Schrodinger 
equation with standard harmonic potential 

K^ +(/32t)2 )^ =a ^' a <^ 2 H)- 

If we set 

= <i>(Hu 

l£N 
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Figure 10. 2D contour plots, in the square [—5, 5] 2 , of an excited 
state solution of the system. In the left figures (two nodal regions) 
we have the 4>\ component corresponding to N\ = N2 = m\ = 
m<z = 1, l\ = 1, I2 = (initial guess), $11 = 50, $22 = 5, $12 = 
(top) and i?i2 = 120 (bottom). In the right figures (no nodal 
regions) we have the corresponding (f>2 component, with l\ = 1% = 
(initial guess). The potentials are centered at the origin. 



where 

fa = (<f>,Hi) L2 = f <Wu 

JR 

the energy functional rewrites as 

tf(*) = 5>#+/ (vw-^^fewu) + W few) 

i£N ^ 1 ' ViSN / 1 jM V/6N / 



and the chemical potential turns into 
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Figure 11. 2D contour plots, in the square [—5, 5] 2 , of an excited 
state solution of the system. In the left figures (four nodal regions) 
we have the <f>\ component corresponding to N-y = N2 = mi = 
m-2 = 1, l\ = I2 = 1 (initial guess), #n = 10, #22 = 5, $12 = 
(top) and $12 = 120 (bottom). In the right figures (no nodal 
regions) we have the corresponding (j>2 component, with l\ = I2 = 
(initial guess). The potentials are centered at the origin. Visibly, 
the supports of the 4>iS segregate around the origin. 



By minimizing E, under the constraint ||0||^2 = N, we look for local minima of 

E(0;A)=£($ + ANV-5>M 
V (eN / 

which solve the system, with k £ N, 

I (a k - a)0 k + jf (v( X ) - u k frn^j +d j v Hk (e ^ n ^j = °> 

\ ten 

We notice that, if is a solution of the above system, then it is immediately seen, 
by multiplying times summing up over k and using (|4.ip . that the Lagrange 
multiplier A equals the chemical potential \i. Next, we truncate to degree L—l and 
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introduce an additional parameter p = 1 in front of the first integral (its usage will 
be clear later), to obtain a corresponding truncated energy functional El{4>', A; p), 
whose local minima solve the system, with < k < L — 1, 

(x K -X)(t> K + pJ^(v(x)-^-^ n k (X> W 'J +d j v Hk =0, 

In order to approximate the integrals, we used a Gauss-Hcrmite quadrature formula 
with 2L — 1 nodes relative to the weight e _2/3 x . Using the tensor basis of the 
Hcrmitc functions, i.e. 

H l (x 1 ,x 2 ) = H^(x 1 )H^(x 2 )e~i^+^ 

the extension to the two-dimensional case is straightforward. In particular, in K 2 , 
7^0,0(^1 j x 2) is the ground eigenstate and r hLi 1 ^ 2 {xi,x 2 ) with any h ^ or l 2 7^ is 
an excited eigenstate of the Schrodingcr equation with standard harmonic potential. 
See Figures [T2l and [T3l representing H.1.0, Jii i, / H 2 .i and 7^2,2- For small coupling 
constants excited states solutions of the GPE system look like these profiles, see 
e.g. Figures [TU1 and [TT1 The extension to a system of any number of equations is not 
difficult, too. In fact, it is sufficient to consider the total energy of the system as the 
functional to be minimized, with a normalization constraint (Lagrange multiplier) 
for each wave function. 



"Hi, 0(^1,2:2) ^14(21,2:2) 




Figure 12. 2D contour plots of Hcrmitc functions with l\ = 1, 
l 2 = 0, Pi = (3 2 = 1 (left picture, one nodal region) and h = l 2 = 1, 
Pi = ft 2 = 1 (right picture, two nodal regions). 



The system is solved by a modified Newton method with backtracking line- 
search, which guarantees global convergence to the ground states. We refer to [6,9] 
and, in particular, to [8] for the details. Here we just mention that only the diagonal 
part of the Jacobian relative to 4>i is computed, thus leading to a dramatic reduction 
of the computational cost for the solution of each linear system. Moreover, the 
initial guess for the Newton iteration is obtained by a continuation technique over p 
and i?, starting from the ground state of the Schrodingcr equation with the standard 
harmonic potential, which corresponds to p = d = 0. The convergence to the 
excited states is not guaranteed, although we observed numerical convergence for 
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Figure 13. 2D contour plots of Hcrmite functions with l-y = 2, 
l 2 = 1, Pi = = 1 (left picture, three nodal regions) and li = 
I2 = 2, /?i = f3 2 = 1 (right picture, four nodal regions). 



the examples reported in the previous section. Of course the case of very large 
values of the coefficients fiij can be treated within the framework of the Thomas- 
Fermi approximation. 
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